Risk Maps

bite <- function(T){
  bite <- (1.67*10^-4) * T * (T- 2.3) * (32.0 - T)^(1/2)
  bite <- ifelse(is.nan(bite), 0, bite)
  return(bite)
}

transmit <- function(T){
  -(2.94*10^-3) * T * (T - 11.3) * (T - 41.9)
}

Figure 2 - Mordecai equation plot

landcover_avg <- readRDS(here::here("risk_maps",
                                   "data", 
                                   "processed_data", 
                                   "landcover_avgs.RData"))

Paired t-tests by land cover type

# Air temperature
t.test(filter(landcover_avg, landcover == "Agriculture")$air_temperature,
       filter(landcover_avg, landcover == "Urban")$air_temperature,
       paired = TRUE)

## 
##  Paired t-test
## 
## data:  filter(landcover_avg, landcover == "Agriculture")$air_temperature and filter(landcover_avg, landcover == "Urban")$air_temperature
## t = -42.019, df = 64, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.084814 -1.895572
## sample estimates:
## mean of the differences 
##               -1.990193

# biting rate
t.test(filter(landcover_avg, landcover == "Agriculture")$biting_rate,
       filter(landcover_avg, landcover == "Urban")$biting_rate,
       paired = TRUE)

## 
##  Paired t-test
## 
## data:  filter(landcover_avg, landcover == "Agriculture")$biting_rate and filter(landcover_avg, landcover == "Urban")$biting_rate
## t = 3.2143, df = 64, p-value = 0.00205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.007554222 0.032364029
## sample estimates:
## mean of the differences 
##              0.01995913

# transmission prob
t.test(filter(landcover_avg, landcover == "Agriculture")$transmission_prob,
       filter(landcover_avg, landcover == "Urban")$transmission_prob,
       paired = TRUE)

## 
##  Paired t-test
## 
## data:  filter(landcover_avg, landcover == "Agriculture")$transmission_prob and filter(landcover_avg, landcover == "Urban")$transmission_prob
## t = -5.0609, df = 64, p-value = 3.767e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.5089083 -0.6548028
## sample estimates:
## mean of the differences 
##               -1.081856

Figure 3- integrated (no temporal variability)

all_pixels_integrate <- readRDS(here::here("risk_maps",
                                         "data", 
                                         "processed_data", 
                                         "all_pixels_integrate.RData"))

coordinates(all_pixels_integrate) <- ~ x + y
gridded(all_pixels_integrate) <- TRUE
air_temp_raster <- raster(all_pixels_integrate, "air_temperature")
biting_raster <- raster(all_pixels_integrate, "biting_rate")
transmission_raster <- raster(all_pixels_integrate, "transmission_prob")

pal <- colorRampPalette(c("darkblue", "blue", "steelblue1", "white", "yellow", "orange", "red"))

plot(air_temp_raster,
       col = pal(50),
       # xlim=c(0,55),
       # ylim=c(-10,50),
       # zlim=limits,
       xlab="Longitude",
       ylab = "Latitude")
title(TeX("Air temperature ($^o$C)"))

plot(biting_raster,
       col = pal(50),
       # xlim=c(0,55),
       # ylim=c(-10,50),
       # zlim=limits,
       xlab="Longitude",
       ylab = "Latitude")
title(TeX("Biting rate ($day^{-1}$)"))

plot(transmission_raster,
       col = pal(50),
       # xlim=c(0,55),
       # ylim=c(-10,50),
       # zlim=limits,
       xlab="Longitude",
       ylab = "Latitude")
title(TeX("Transmission probability (%)"))

Figure 3 - maps for each time of day

Air Temperature

Figure S3 - maps of errors

Maps of errors integrated

get_error <- function(temp_raster, risk_function){

  max_error_by_pixel <- function(temperature){
    if (is.na(temperature)){
      return(NA)
    } else {
      lower_bound <- temperature - (STANDARD_ERROR/(65-1))*1.96
      upper_bound <- temperature + (STANDARD_ERROR/(65-1))*1.96
  
      max <- optimize(risk_function, interval=c(lower_bound, upper_bound), maximum=TRUE)$objective
      min <- optimize(risk_function, interval=c(lower_bound, upper_bound), maximum=FALSE)$objective
  
      return(max - min)
    }
  }

  values <- raster::values(temp_raster)
  errors <- sapply(values, max_error_by_pixel)
  values(temp_raster) <- errors
  return(temp_raster)
}

plot(get_error(air_temp_raster, bite),
       col = pal(50),
       # xlim=c(0,55),
       # ylim=c(-10,50),
       # zlim=limits,
       xlab="Longitude",
       ylab = "Latitude")
title(TeX("Biting rate uncertainty ($day^{-1}$)"))

plot(get_error(air_temp_raster, transmit),
       col = pal(50),
       # xlim=c(0,55),
       # ylim=c(-10,50),
       # zlim=limits,
       xlab="Longitude",
       ylab = "Latitude")
title(TeX("Transmission probability uncertainty (%)"))

Biting rates errors

## NULL

transmission rates errors – CODE WRONG. NEEDS CHECK.

Figure 4 - Distribution of temporal averages over crop types

Integrated values over whole day for pixels close to each other

Because each ECOSTRESS image has a slightly different grid, it is difficult to match the pixels to each other. We do this by resampling all images to the first image.

pixels <- readRDS(here::here("risk_maps",
                                         "data", 
                                         "processed_data", 
                                         "all_pixels_integrate_w_landcover.RData")) %>% 
  filter(!is.na(air_temperature))

ag_types <- c("Vegetable", "Orchard", "Fruit", "Field Crop", "Uncultivated", "Urban")
pixels$landcover <- factor(pixels$landcover, levels = c("Urban",  "Uncultivated", "Vegetable", "Fruit", "Orchard", "Field Crop", "Other"))

filter(pixels, landcover %in% ag_types) %>%
  ggplot(aes(x = air_temperature, col= landcover, alpha = 0.2)) +
  geom_density() +
    guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
      xlab(TeX("Air temperature ($^o$C)")) +
      ylab("Density") +
  xlim(22.5, 29.2) +
  # labs(title = "Average air temperature distribution") +
  theme_bw() +
  theme(legend.position = c(0.15, 0.7), legend.title = element_blank()) 

filter(pixels, landcover %in% ag_types) %>%
  ggplot(aes(x = transmission_prob, col= landcover, alpha = 0.2)) +
  geom_density() +
    guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
      xlab("Transmission probability (%)") +
      ylab("Density") +
  xlim(13, 17.3) +
  # labs(title = "Average transmission probability distribution") +
  theme_bw() + 
  guides(color = FALSE)

filter(pixels, landcover %in% ag_types) %>%
  ggplot(aes(x = biting_rate, col= landcover, alpha = 0.2)) +
  geom_density() +
    guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
      xlab(TeX("Biting rate ($day^{-1}$)")) +
      ylab("Density") +
  xlim(.11, .225) +
  # labs(title = "Average biting rate distribution") +
  theme_bw() + 
  guides(color = FALSE)

# original submission: Instead of just looking at a single daytime frame's pixels, look at all measurements from 11 to 16h (15 frames)

# day_pixels <- readRDS(here::here("risk_maps",
#                                "data", 
#                                "processed_data", 
#                                "day_pixels.RData"))
# 
# ag_types <- c("Vegetable", "Orchard", "Fruit", "Field Crop", "Uncultivated")
# day_pixels$landcover <- factor(day_pixels$landcover, levels = c("Urban",  "Uncultivated", "Vegetable", "Fruit", "Orchard", "Field Crop", "Other"))
# 
# filter(day_pixels, landcover %in% ag_types) %>%
#   ggplot(aes(x = air_temperature, col= landcover, alpha = 0.2)) + 
#   geom_density() + 
#     guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
#       xlab(TeX("Air temperature ($^o$C)")) + 
#       ylab("Density") + 
#   xlim(25, 36.7) + 
#   labs(title = "Daytime air temperature distribution") + 
#   theme_bw()
# 
# filter(day_pixels, landcover %in% ag_types) %>%
#   ggplot(aes(x = transmission_prob, col= landcover, alpha = 0.2)) + 
#   geom_density() + 
#     guides(alpha = FALSE, col = guide_legend(title = "Crop category")) +
#       xlab("Transmission probability (%)") + 
#       ylab("Density") + 
#   xlim(14.5, 20) + 
#   labs(title = "Daytime transmission probability distribution") + 
#   theme_bw()

f-test on different crop types

f_air_temp <- aov(air_temperature ~ landcover, data = filter(pixels, landcover %in% ag_types))
f_trans_prob <- aov(transmission_prob ~ landcover, data = filter(pixels, landcover %in% ag_types))

summary(f_air_temp)

##                 Df Sum Sq Mean Sq F value Pr(>F)    
## landcover        5 333260   66652  120990 <2e-16 ***
## Residuals   448948 247321       1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(f_trans_prob)

##                 Df Sum Sq Mean Sq F value Pr(>F)    
## landcover        5  87220   17444  145652 <2e-16 ***
## Residuals   448948  53768       0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Figure 5 - urban and ag diurnal time series

Jensen’s inequality test

####(space averages)

All <-filter(landcover_avg, landcover == "All")

bite <- function(T){
  bite <- (1.67*10^-4) * T * (T- 2.3) * (32.0 - T)^(1/2)
  bite[is.nan(bite)] <- 0
  return(bite)
}

transmit <- function(T){
  -(2.94*10^-3) * T * (T - 11.3) * (T - 41.9)
}

All$avg_bite <- All$air_temperature %>% bite
All$avg_transmit <- All$air_temperature %>% transmit

t.test(x = All$biting_rate, y = All$avg_bite, paired = TRUE)

## 
##  Paired t-test
## 
## data:  All$biting_rate and All$avg_bite
## t = 2.5515, df = 64, p-value = 0.01313
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.001369638 0.011251819
## sample estimates:
## mean of the differences 
##             0.006310729

t.test(x = All$transmission_prob, y = All$avg_transmit, paired = TRUE)

## 
##  Paired t-test
## 
## data:  All$transmission_prob and All$avg_transmit
## t = -7.7387, df = 64, p-value = 9.314e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3090675 -0.1822379
## sample estimates:
## mean of the differences 
##              -0.2456527

####(time averages)

all_pixels_integrate <- readRDS(here::here("risk_maps",
                               "data", 
                               "processed_data", 
                               "all_pixels_integrate.RData")) %>% filter(!is.na(air_temperature)) #there are 3751 pixels on the edge that do not have values for every image, and these we remove. 

all_pixels_integrate$avg_bite <- all_pixels_integrate$air_temperature %>% bite 
all_pixels_integrate$avg_transmit <- all_pixels_integrate$air_temperature %>% transmit

t.test(x = all_pixels_integrate$biting_rate, y = all_pixels_integrate$avg_bite, paired = TRUE)

## 
##  Paired t-test
## 
## data:  all_pixels_integrate$biting_rate and all_pixels_integrate$avg_bite
## t = -4384.6, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.09894369 -0.09885528
## sample estimates:
## mean of the differences 
##             -0.09889949

t.test(x = all_pixels_integrate$transmission_prob, y = all_pixels_integrate$avg_transmit, paired = TRUE)

## 
##  Paired t-test
## 
## data:  all_pixels_integrate$transmission_prob and all_pixels_integrate$avg_transmit
## t = -3129.6, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.069540 -3.065697
## sample estimates:
## mean of the differences 
##               -3.067618

role of both spatial and temporal resolution

all_pixels_integrate <- readRDS(here::here("risk_maps",
                               "data", 
                               "processed_data", 
                               "all_pixels_integrate.RData")) %>% filter(!is.na(air_temperature)) #there are 3751 pixels on the edge that do not have values for every image, and these we remove. 

avg_bite <- all_pixels_integrate$air_temperature %>% mean(na.rm = TRUE) %>% bite 
avg_transmit <- all_pixels_integrate$air_temperature %>% mean(na.rm = TRUE) %>% transmit

t.test(x = all_pixels_integrate$biting_rate, mu = avg_bite)

## 
##  One Sample t-test
## 
## data:  all_pixels_integrate$biting_rate
## t = -4372.7, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0.2512659
## 95 percent confidence interval:
##  0.1492290 0.1493204
## sample estimates:
## mean of x 
## 0.1492747

t.test(x = all_pixels_integrate$transmission_prob, mu = avg_transmit)

## 
##  One Sample t-test
## 
## data:  all_pixels_integrate$transmission_prob
## t = -5128.9, df = 964748, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 18.21251
## 95 percent confidence interval:
##  15.04580 15.04822
## sample estimates:
## mean of x 
##  15.04701

bitemax <- optimize(bite, interval=c(0, 40), maximum=TRUE)$objective
transmax <- optimize(transmit, interval=c(0, 40), maximum=TRUE)$objective

Plot Jensen’s inequality for spatial stuff through time: not included

ggplot(All) + geom_point(aes(x = avg_bite - biting_rate, y = hour(dt) + minute(dt)/60, color = avg_bite))

ggplot(All) + geom_point(aes(x = avg_transmit - transmission_prob, y = hour(dt) + minute(dt)/60, color = avg_transmit))

Figure Sx: Distribution of ECOSTRESS images included in the study by month and hour of day.

lstfiles <- list.files(path = here::here("risk_maps", 
                                         "data", 
                                         "raw_data", 
                                         "ECOSTRESS",
                                         "all_ims",
                                         "LST"), pattern = "tif", recursive = TRUE, full.names = FALSE)

getdate <- function(file){
  
  nums <- str_extract(file, regex('(?<=doy)[0-9]+'))
  year <- substr(nums, 1, 4)
  doy <- substr(nums, 5, 7) %>% as.numeric()
  hhmmss <- substr(nums, 8, 13)
  date <- as.Date(doy - 1, origin = paste0(year, "-01-01"))
  dt <- ymd_hms(paste(date, hhmmss), tz = "UTC") %>% with_tz("America/Los_Angeles")
  
  return(dt)
}

all_dates <- do.call(c, lapply(lstfiles, getdate))

files <- list.files(path = here::here("risk_maps", 
                                      "data",
                                      "raw_data",
                                      "ECOSTRESS",
                                      "filtered_ims"), 
                    pattern = "tif", 
                    recursive = TRUE, 
                    full.names = FALSE)

getfiltdate <- function(file){
  
  date <- ymd(substr(file, 14, 23))
  hhmmss <- str_extract(file, regex('[0-9]{2}:{1}[0-9]{2}:{1}[0-9]{2}'))
  dt <- ymd_hms(paste(date, hhmmss), tz = "America/Los_Angeles")
  
  return(dt)
}

filtered_dates <- do.call(c, lapply(files, getfiltdate))

data.frame(dt = all_dates, included = ifelse(all_dates %in% filtered_dates, "Included", "Not included")) %>%
  ggplot() +
  geom_point(aes(x = month(dt) + day(dt)/30 + hour(dt)/(24*30), y = hour(dt) + minute(dt)/60, shape = as.factor(year(dt)), color = included)) +
  scale_color_manual(values = c("red", "lightblue2")) + 
  xlab("Month") + 
  ylab("Hour")  +
  theme_bw() +
  theme(legend.title = element_blank())